Violation of hyperbolicity in a diffusive medium with local hyperbolic attractor 
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Departing from a system of two non-autonomous amplitude equations, demonstrating hyperbolic 
chaotic dynamics, we construct a ID medium as ensemble of such local elements introducing spatial 
coupling via diffusion. When the length of the medium is small, all spatial cells oscillate syn- 
chronously, reproducing the local hyperbolic dynamics. This regime is characterized by a single 
positive Lyapunov exponent. The hyperbolicity survives when the system gets larger in length so 
that the second Lyapunov exponent passes zero, and the oscillations become inhomogeneous in space. 
However, at a point where the third Lyapunov exponent becomes positive, some bifurcation occurs 
that results in violation of the hyperbolicity due to the emergence of one- dimensional intersections 
of contracting and expanding tangent subspaces along trajectories on the attractor. Further growth 
of the length results in two-dimensional intersections of expanding and contracting subspaces that 
we classify as a stronger type of the violation. Beyond of the point of the hyperbolicity loss, the 
system demonstrates an extensive spatiotemporal chaos typical for extended chaotic systems: when 
the length of the system increases the Kaplan- Yorke dimension, the number of positive Lyapunov 
exponents, and the upper estimate for Kolmogorov-Sinai entropy grow linearly, while the Lyapunov 
spectrum tends to a limiting curve. 
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Introduction 

One of the central concepts in mathematical theory of 
dynamical systems relates to hyperbolic strange attrac- 
tors. Tangent space of each point of such an attractor 
splits into expanding and contracting subspaces, and this 
splitting is invariant. Dynamics on a hyperbolic attrac- 
tor is structurally stable, i.c, is insensible to variations 
of parameters. It manifests strong stochastic properties 
and allows detailed theoretical analysis P, Q ■ 

During the last 40 years hyperbolic attractors were 
considered rather as idealized model of perfect chaos. 
Though some artificial systems with hyperbolic attrac- 
tors were known, they were useless for practical appli- 
cations because of complicated construction. Recently, 
a realistic system was suggested and implemented as 
electronic device, dynamics of which in stroboscopic de- 
scription is associated with attractor of Smalc- Williams 
type d, B- Attractor of this system is hyperbolic as 
proven numerically by verification of the cone crite- 
rion Q . In paper [y] the amplitude equation for this sys- 
tem was studied, and the hyperbolicity was also proven 
by the method of cones. (Some other models based on 
the same principle are considered in Refs. @, 11, & EH ) 

Traditionally, studies of hyperbolic dynamics are 
mostly concentrated on low dimensional systems. Many 
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topics concerning spatiotemporal chaos, though at- 
tracted a lot of interest, remain open [111 ]. In this pa- 
per we address a problem of survival of hyperbolicity of 
a spatiotemporal system when the length of the system 
grows. We consider a ID extended system composed of 
local elements possessing a hyperbolic attractor that is 
based on the amplitude equations from The spatial 
coupling is introduced via diffusion. In fact, a system we 
study is a set of two coupled non-autonomous Ginzburg- 
Landau equations of special form. 

We are aware of two numerical methods for reliable 
verification of hyperbolicity. The first one is the method 
based on the cone criterion ||, which employs directly 
the rigorous theorem, and, hence, looks preferable. Un- 
fortunately, the method is appropriate only for low- 
dimensional systems, while its extension to systems of 
many degrees of freedom seems to be abundantly sophis- 
ticated. The second method is based on a recently sug- 
gested routine of computing of covariant Lyapunov vec- 
tors [12]. These vectors are associated with Lyapunov 
exponents and indicate directions of contracting and ex- 
panding manifolds at each point of an attractor. If these 
vectors are known, angles between each contracting and 
each expanding direction can be computed and the min- 
imal one can be found. The attractor is interpreted as 
non- hyperbolic, if the distribution of these angles does 
not vanish at the origin. In fact, this is only a suffi- 
cient condition because the converse is not true. In the 
present paper we apply more subtle approach based on 
computation of so called principal angles fl3| . that al- 
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lows to detect a tangency of two arbitrary vectors from 
contracting and expanding subspaces. 

The paper is organized as follows. In Sec. U we intro- 
duce the system and briefly discuss its local dynamics. 
Also we describe a numerical method applied to find so- 
lutions to the system. Sec. [TT| represents linear stability 
analysis. The critical length of the system is determined 
where a spatially homogeneous solution becomes unsta- 
ble with respect to non-uniform perturbation. Sec. IIIII is 
devoted to illustrations of spatiotemporal dynamics. The 
main part of the paper is Sec. IIV1 where we develop the 
Lyapunov analysis. We discuss distributions of minimal 
angles between contracting and expanding tangent sub- 
spaces on the attractor. Also, dependencies of Lyapunov 
exponents, Kaplan- Yorke dimension and Kolmogorov- 
Sinai entropy on the length of the system are considered. 
In Sec. |V] we summarize the obtained results and outline 
perspectives for further investigations. 

I. THE MODEL AND NUMERICAL METHOD 

Let us start with a physical model, demonstrating hy- 
perbolic dynamics, suggested by Kuznetsov in Ref. [3|. 
The model consists of two coupled non-autonomous van 
der Pol oscillators that are parametrically influenced by 
an external periodic force. The oscillators become active 
turn by turn, and pass the excitation each other in such 
way that the phase of oscillations is doubled after each 
period of the forcing. In Ref. 6] the amplitude equations 
for this system are derived that read 

a = Acos(2nt/T)a - \a\ 2 a - ieb, 
b = -Acos(2nt/T)b- \b\ 2 b~iea 2 . '"' ' 

In this paper we study a spatially extended analog of 
these equations, supplying them with the second spatial 
derivatives. 

So, we consider two coupled non-autonomous 
Ginzburg-Landau equations: 

d t a = A cos(27rf/T)a - \a\ 2 a -ieb + d£a, 
d t b = -Acos(2nt/T)b - \b\ 2 b - ie a 2 + d 2 x b. 

Here a = a{x, t) and 6 = b[x, t) are complex dynami- 
cal variables whose behavior is the subject of interest. 
Coefficients at linear terms undergo periodic variation 
with period T and amplitude A. The parameter modu- 
lation takes place in opposite phase for a and b. When 
the first subsystem is excited, the second one is relaxed 
and vice versa. The forcing is supposed to be slow, i.e., 
the half period T/2 is much longer then a transient time 
of the excitation. The second terms in the right-hand 
parts of the equations provide saturation of instabilities 
in the excited subsystems. Additionally, there are terms 
responsible for coupling between the components a and 
b] the intensity of the coupling is controlled by e. The 
coupling is asymmetric, being quadratic from a to & and 



linear in the inverse direction. Finally, the last terms in 
the right-hand parts introduce diffusion, that is respon- 
sible for the spatial distribution of local oscillations. The 
diffusion coefficients of the subsystems are equal to 1. We 
study the system in a limited spatial domain < x < L. 
The boundary conditions are 

( 9 * a )\x*4},L = ( d * h )\x=0X = °- ( 3 ) 

Let us briefly discuss a local dynamics of the system. 
Consider Eqs. |T]). (A more detailed study can be found 
in Ref. (&].) Due to the presence of periodic forcing in 
([T]), it is natural to introduce a stroboscopic map: we 
split the continuous time into steps of length T, and con- 
sider a sequence of states of the system at the beginnings 
of these steps. Define phases within the interval [0,2-71"): 
4> = arga, ip = arg&. Suppose at some instant the first 
oscillator is excited, and its amplitude \a\ is high. Then, 
the second one is suppressed, and its amplitude |6| is 
small. The coefficients in (JXJ) are real, except the cou- 
pling term. It means that the phases can vary only as 
a result of interaction between subsystems. But, when 
a is excited, |6| is small, and its action on a is negligi- 
ble. Thus, the phase of a remains approximately con- 
stant during the excitation stage. On the contrary, the 
influence of the excited a on the suppressed b is strong. 
The coupling term is proportional to a 2 . It means that 
after the half period T/2 at the threshold of its own exci- 
tation the oscillator b inherits a doubled phase of a (also 
the phase gets a shift — n/2 because of the imaginary unit 
at the coupling term). Now the roles of the subsystems 
are exchanged. The phase of b remains constant when 
this subsystem is excited and at the end, after the other 
T/2, the phase is returned back to a through a linear cou- 
pling term (also with the shift —ir/2). As a result, the 
first oscillator a doubles its phase during the period T. 
This discussion allows to write down a map for a series of 
phases <j> n = arga(nT) that are measured over the time 
step T: 

<f>n+i = 2(j) n - n mod 2n. (4) 

Up to a constant term (that can be eliminated by a 
shift of the origin of the phase) this map coincides with 
the well known Bernoulli map 

El El- 

It demonstrates 

chaotic dynamics, and the chaos is homogeneous: a rate 
of exponential divergence of two close trajectories is iden- 
tical at each point of the phase space, being equal to In 2. 

Getting back to the continues system ([T]), we estimate 
its largest Lyapunov exponent as 

A = ln2/T. (5) 

The described mechanism of phase doubling presumes a 
hyperbolic nature of the dynamics of (Q} . The numerical 
verification of the cone criterion, that has been preformed 
in [6j, confirms this. 

Before starting an analysis of the system @, let us 
discuss a numerical method applied to find its solutions. 
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Formally, our equations can be classified as parabolic 
PDE. Typical recommendation of handbooks for such 
equations is the Crank-Nicolson method which is abso- 
lutely stable and provides the second order of local ap- 
proximation both in space and in time. This method is 
semi- implicit, i.e., a solution at a new level tk+i is ex- 
pressed via previous solution at tk as a set of algebraic 
equations, so that values from all spatial points on both 
levels are involved into this equation set. If PDE is linear, 
these equations are linear too. But application of this ap- 
proach to non-linear PDEs, like ours, gives rise to a set of 
non-linear algebraic equations that requires much more 
computational efforts. Usually, one simplifies the prob- 
lem by neglecting terms, being non-linear with respect 
to unknown variables. The resulting numerical scheme is 
semi-implicit for linear part of initial PDE and explicit 
for non-linear part. Unfortunately, this simplified "quasi 
Crank-Nicolson" method is not absolutely stable. Some- 
times everything goes fine, but sometimes, usually when 
the system is far beyond the instability threshold, the 
solution diverges. In this paper we do not neglect the 
non-linearity and develop a true semi-implicit scheme. 
At each time step we solve a set of non-linear equations 
via the Newton-Raphson iterations. The seed for the it- 
erations is found from the mentioned simplified method. 
The iterations converge very fast. Normally, it takes 2 or 
3 repetitions to solve the non-linear equations with the 
accuracy 10 -5 or even better. The idea of the described 
method can be found in books on numerical analysis, e.g., 
[HI, G3]- Though the method is a bit complicated, this is 
compensated by its high accuracy and stability. 

Below different characteristic values are calculated as 
functions of the length of the system L. Varying L, we 
need to choose some strategy of simultaneous variation 
of parameters of a numerical mesh. One way is to keep 
constant number of points of the mesh N and compute 
space step as Ax = L/(N — 1). The other way is to fix 
the step Ax and find N for each L as N = 1 + \L/Ax] , 
where [•] means ceiling (to get a consistent numerical 
scheme one also needs to adjust actual value of Ax for 
the equality Ax(N — 1) = L to fulfill). In our simula- 
tions we always keep constant N. This strategy seems to 
be preferable because the number of degrees of freedom 
of the numerical model remains constant; obviously, it 
is equal to 2N (traditionally defined as a half of a total 
order of the set of differential equations). So, we can be 
sure that phenomena, observed when L is varied, emerge 
due to a transformation of an inner structure of attrac- 
tor, and they can not be attributed to just an extensive 
increase of degrees of freedom. The time step At can 
be either constant or attached to Ax. When At is suf- 
ficiently small, these two ways produce identical results. 
We shall hold the time step at At w 0.01. (Addition- 
ally, a small adjustment is also made to fit an integer 
number of steps into the observation interval). Though 
this is redundantly small value to obtain solutions to the 
system ©, but this is needed to estimate correctly its 
minor Lyapunov exponents. 



II. LINEAR STABILITY ANALYSIS 

Standard linear stability analysis of autonomous spa- 
tially extended active system requires a consideration of 
small perturbations to a homogeneous steady state. Ex- 
istence of perturbation modes with positive growth rates 
indicates the instability of the homogeneous state. Our 
system does not have a steady state, and its dynamics is 
chaotic in time. Oscillations can be either homogeneous 
or irregular in space. Our aim is to find the conditions for 
a transition from one regime to another, utilizing ideas 
of the standard analysis. 

Suppose that the system is infinite in space and its ini- 
tial state is uniform. Prepared in this way, the system 
obviously demonstrates homogeneous oscillations; at any 
spatial point the dynamics can be described by the ODE 
system ([1]). Let us consider an inhomogeneous pertur- 
bation to these oscillations. We need to seek a solution 
composed as a sum of a homogeneous part, say, ao(i) 
and bo(t), and a sinusoidal mode of perturbation with 
real wave number k and real growth rate a(k). The sys- 
tem is chaotic; thus, instead of usual assumption of time 
periodicity of small perturbation, we introduce small am- 
plitudes a(t) and b(t) and require them neither grow nor 
decay, in average. It means that there exist two con- 
stants, < K < M < oo, such that K < \a(t)\ < M and 
K < \b(t)\ < M for t > 0. So, we set: 

a(x,t) = a Q (t) + a(t)e^ t - ik3: , 

b(x,t) = b (t) + b(t)e a W t - ikx . ((>1 

After substitution ([6]) to ([2]), we exclude non-linear terms 
in a and b, supposed to be small, and obtain a set of linear 
ODE for complex amplitudes of perturbation: 

a = (A cos(27rf/T) - A )a - 2|a | 2 a - a„a* - ieb, 

b = {-A cos(27rf/T) - \ )b - 2\b a \ 2 b - bib* - 2iea a, 

(7) 

where Ao = k 2 + <r(k) and asterisk denote the complex 
conjugation. A value of Ao controls growth or decay of a 
solution. Because a{t) and b(t) should be bounded at any 
k, A does not depend on k. One can easily check that ([7]) 
also describes small perturbation to an orbit of {T]). It 
means that the conditions on a(t) and b(t) are fulfilled 
when Ao is equal to the largest Lyapunov exponent of |T|) . 
Thus, we can write 

a(fc) = A - k 2 . (8) 

Relation (JSJ) can be verified by direct numerical com- 
putations of a{k). For this purpose we substitute Ao — > 
cr(fc) + k 2 to ||7J and set there a(k) = 0. It meas that 
now the amplitudes a and b are allowed to grow or decay, 
so that the rate will be equal to a(k). Given k, we find 
<r(k) employing the algorithm of computing of the largest 
Lyapunov exponent [18| . System ([7]) is initialized with 
a unit vector, and then solved together with ([T]) on one 
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Figure 1: Growth rate a(k) of non-uniform perturbation to a 
homogeneous solution of the system ([2]) . The solid curves are 
the graphs of ©: the upper one is (ln2)/5 — k 2 , the lower 
one is <r(0) — k 2 , where <j(0) = 0.00842. Crosses, circles and 
stars represent numerical data, computed at parameters that 
are shown in the legend. 

period T. After that, a norm of the vector-solution of 
© is found and stored, and the vector itself is normal- 
ized. When this procedure is repeated for a sufficiently 
long time, the averaged logarithms of the collected norms 
determine the sought a(k). The results are shown in 
Fig. [TJ Solid lines represent theoretical er(fc) ©. The 
upper one corresponds to a hyperbolic chaos in {T]) and 
Ao is found according to ([5]) ■ The lower curve also corre- 
sponds to chaotic oscillations of ©, that are, however, 
non- hyperbolic. In this case we substitute a computed 
value of <t(0) to © instead of Ao- Numerical data fit well 
the theoretical curves. As follows from © and ©, cr(fc) 
does not depend on the A in the regime of hyperbolic 
chaos. Numerical verification confirms this. 

Linear modes described by © are influenced para- 
metrically by a chaotic force. It means that all modes 
with positive cr(fc) can grow simultaneously giving rise 
spatiotemporal chaos. The spectrum of lineally unstable 
modes with a(k) > can be found from ©. These modes 
lay within the interval of wave numbers < k < k\i n , 
where 

fclin = \/Ao. (9) 

If the system © is bounded by the length L, the spec- 
trum of modes allowed by the boundary conditions © 
is 

k n = nn/L, n= 1,2,3,.... (10) 

When L is small, so that k\ > ku n , there are no unstable 
eigenmodes and the system demonstrates homogeneous 
oscillations. Spatial structure emerges above the critical 
point which can be found from the condition k\ = k\i n : 
L c = 7r/vAo- Below we put attention to the case when 
the local dynamics is hyperbolic. The Lyapunov expo- 
nent Ao in this case is given by © and the critical length 
reads: 

L c = tt^T/ In 2. (11) 



Figure 2: Spatiotemporal dynamics of ([2]) at L = 8. Grey 
levels indicate values of Sfta. A — 3, T = 5, e = 0.05. 



III. SPATIOTEMPORAL DYNAMICS 

Let us consider some illustrations of spatiotemporal 
dynamics of the system © . Figure [5] represents homo- 
geneous oscillations. In this and subsequent figures the 
space coordinate is horizontal, time is directed vertically 
and grey levels indicate values of 3?a as shown by gra- 
dient bars at the right edges of the diagrams. Layers 
9?a(x) are plotted at successive steps t n — nT . Critical 
length, according to (JTTJ) , is L c sa 8.44. The length of 
the system in Fig. [2] is less then the critical value, L = 8. 
Hence, after a short transient time, it settles in a regime 
of homogeneous oscillations. 

In Fig. [3] the length L = 10 is larger than L c . The 
first eigenmode cos(fcij;) falls into the instability domain 
and grows, destroying the homogeneity. The first mode 
contains one half of the period of cosine, so if a maximum 
is at the left edge of the system, a minimum appears 
at the right edge and vice versa. Careful inspection of 
Fig. [3] confirms this conclusion. If a horizontal stripe, 
representing $ia(x) at a certain time step, is white at the 
left edge, it becomes dark at the right edge. 

The result of further increase of the length up to 
L = 500 is shown in Fig. As here a lot of eigenmodes 
satisfy the condition k n < kn n , they are exited and pro- 
duce a rich and complicated structure. It is interesting to 
note that it reminds a structure generated by a cellular 
automata of Wolfram's class 3 fl9ll. 



IV. LYAPUNOV ANALYSIS 

Lyapunov exponents are average rates of expansion or 
contraction in the tangent space on an attractor. They 
characterize the sensitivity of motion to small perturba- 
tions; an attractor with a positive exponent is chaotic. 
Also, it is important to know a mutual orientation of ex- 
panding and contracting directions in the tangent space 
at each point of the attractor. This information can be 
provided by covariant Lyapunov vectors (T3 |. If there is 
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Figure 3: Same as Fig. [3] but at L = 10. 




Figure 4: L = 500. 



a well defined split of the tangent space into contracting 
and expanding subspaces, the dynamics is hyperbolic. 
On the contrary, the dynamics is non-hyperbolic when 
couples of collinear vectors from contracting and expand- 
ing subspaces can be encountered with a non-zero prob- 
ability. 

In these section we compute covariant Lyapunov vec- 
tors and perform a verification of hyperbolicity of the at- 
tractor of ([2]). Also we analyze Lyapunov exponents for 
the system ([2]) as well as related to them Kaplan- Yorkc 
dimension and Kolmogorov- Sinai entropy. 



A. Verification of hyperbolicity at different lengths 
of the system 

To verify the hyperbolicity one needs to analyze ex- 
panding and contracting directions in the tangent space 
on an attractor. These directions can be found in a form 
of covariant Lyapunov vectors [l^] • The method of com- 
putation of these vectors is briefly described in Appendix. 

When the covariant Lyapunov vectors are computed 
at some point of the attractor, the simplest way to verify 
the hyperbolicity is to estimate angles between each cou- 
ple of expanding and contracting vectors and find the 
smallest one. Collecting the smallest angles for suffi- 



ciently many points, one obtains a sufficient condition 
for non- hyperbolicity: the attractor is non-hyperbolic if 
zero angle can be encountered with non-zero probability. 
But the converse is not true. The covariant Lyapunov 
vectors may not be collinear themselves, but the loss of 
hyperbolicity still can take place due to a tangency of 
some other couple of vectors from contracting and ex- 
panding subspaces. To take this situation into account, 
a more subtle approach should be used. 

Let us suppose that at some point of the attractor we 
have n s covariant Lyapunov vectors spanning the con- 
tracting tangent subspace iS and n u vectors that span 
the expanding subspace hi. It is natural to assume that 
n s > n u . Consider unit vectors s € S and u € hi and 
find among them a couple S\ and U\ that produces the 
largest inner product. Arc cosine of sjui is the smallest 
angle between subspaces, that is denoted as 6±. Then we 
seek for unit vectors S2 and U2 that again produce the 
largest inner product but with additional requirement to 
be orthogonal to s% and u\, respectively. Arc cosine of 
their inner product is denoted as 02- Proceeding with 
this procedure, we obtain n u angles, 



< 6>i < . . . < 6»„ u < tt/2, 



(12) 



that are called the principal angles. Corresponding vec- 
tors Si and Ui are called the principal vectors. The formal 
definition of the principal angles and vectors is the fol- 
lowing 



where 



cos Ok = max max s u 

s£S u&A 



s T s — u T u = 1, 
s T Si = 0, u T Ui = 0, i = 1, 



,fc- 1. 



(13) 



(14) 



The algorithm of computation of the principal angles is 
discussed in Appendix. 

Vanish of the principal angles indicate a tangency be- 
tween contracting and expanding subspaces and violation 
of the hyperbolicity. A necessary and sufficient condition 
for the loss of hyperbolicity is appearance of such distri- 
bution of B\ on the attractor that it has a non-zero value 
at the origin. If a system has many degrees of freedom, 
several smallest principal angles can vanish simultane- 
ously, that means that several couples of contracting and 
expanding vectors merge. A number of such angles de- 
fines the dimension of the tangency. A necessary and 
sufficient condition for the n-dimensional tangencies is a 
non-zero probability of vanish of the sum of first n prin- 
cipal angles. 

Figure [5] represents distributions of 9i for the sys- 
tem ^ at different lengths L. The equations have been 
solved at At » 0.01 and Ax = L/(N - 1), where N 
is the number of points of a numerical mesh. N = 51 
for all L, except L = 60 where N = 101. The distri- 
butions have been computed with resolution 300 points. 
For each distribution 10 trajectories of the length 600T 
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Figure 5: Distributions of minimal angles Q\ between con- 
tracting and expanding tangent subspaces of the attractor 
of (J2]). A = 3, T = 5, e = 0.05. The logarithmic scale is 
used along the ordinate axis. The dashed lines in the panels 
L — 17, and L = 30 are obtained via least squares fit. There 
is one positive Lyapunov exponents at L = 8, two at L = 10 
and L = 15, three at L — 17, five at L = 30, and eleven at 
L = 60. Observe the violation of hyperbolicity at L > 15, i.e, 
when the third Lyapunov exponent becomes positive. 



have been processed with the interval T/30 between re- 
normalizations and orthogonalizations (see Appendix for 
details). In the course of the backward iterations, a time 
interval 500T is omitted as transient, and then the an- 
gles are computed on the interval 100T. Thus, totaly 
3000 angles for each trajectory have been stored. The 

distributions have been normalized, J? P{9\) = 1. 

The upper curve L = 8 in Fig. [5] corresponds to a 
spatially homogeneous case. The angles are very well lo- 
calized. Thus, the hyperbolic dynamics is observed that 
corresponds to the hyperbolic dynamics of the ODE sys- 
tem (JTJ) - The second curve L = 10 represents the case of a 
weak inhomogeneity, when the system is not far above the 
critical point L c . Observe that the distribution becomes 
much more smooth, compared to the homogeneous case. 
It means that different configurations of contracting and 
expanding subspaces are encountered with almost equal 
probabilities. The distribution is still separated well from 
the origin, i.e., the attractor remans hyperbolic. This is 
also the case for the next distribution at L = 15. This 
distribution is even more flat than the previous one, and 
also it is separated well from the origin. Notice that there 



are two positive Lyapunov exponents both at L — 10 and 
at L = 15. The picture becomes dramatically different at 
L = 17, when the third Lyapunov exponent becomes pos- 
itive. The distribution occupies almost the whole range of 
angles and has non-zero value at origin. The former indi- 
cates that the attractor becomes non-hyperbolic. More- 
over notice that in the logarithmic scale the curve decays 
linearly from the origin. It means that the most part of 
the distribution is described by an exponential function. 
Similar behavior is observed at L = 30: the most part 
of the curve obeys the exponential law. The exponents, 
that are equal to the slopes of the dashed approximating 
lines, are —0.72 at L = 17, and —1.81 at L = 30, i.e., 
their absolute values grow with L. When L gets larger, 
as in the panel for L = 60, the distribution undergoes a 
transformation. It acquires an extended sloping segment 
near the origin, while the other part of the distribution 
becomes more or less flat, on average. The attractor 
remains non- hyperbolic, and, moreover, the probability 
to encounter the tangency of contracting and expanding 
subspaces becomes larger. 

We can assume that the reorganization of the structure 
of distribution, that occurs between L — 30 and L = 60, 
is associated with emergence of the two-dimensional tan- 
gencies of contracting and expanding subspaces. Figure[|5] 
demonstrates distributions of two first principal angles 
(6 1 + #2)/2. The curve at L — 17 is separated well from 
the origin, so that no two-dimensional tangencies take 
place. At L = 30 the curve approaches zero much closer. 
Finally, the curve at L = 60 touches the ordinate axis 
confirming the presence of the two-dimensional tangen- 
cies. 

Figure [7| reproduces the observed scenario at some 
other set of parameters. In panel (a) we can see that 
the attractor is hyperbolic with two positive Lyapunov 
exponents, curve L — 10, while emergence of the third 
one results in the violation of the hyperbolicity, curve 
L = 11. Similar to the case presented in Fig. the 
distribution right above the violation point is basically 
exponential, curve L = 11, while the further growth of 
L results in the transformation of the distribution, curve 
L = 60. Figure [7(b) indicates the emergence of the two- 
dimensional tangencies in this case: the distributions of 
(6 1 + 6*2 )/2 approach the origin as L grows and touch it 
at L = 60. 

So, we observe that the growth of L first results in the 
violation of hyperbolicity due to one-dimensional tangen- 
cies of contracting and expanding subspaces, and then 
gives rise to two-dimensional tangencies between these 
subspaces. It is natural to suggest, that the tangencies 
of higher dimensions also arise at appropriate lengthes 
of the system. The violation of hyperbolicity is accom- 
panied by the emergence of the third positive Lyapunov 
exponent. Let us denote the point where the third Lya- 
punov exponent passes zero as Li- We suspect that the 
loss of hyperbolicity takes place exactly at L = L2, and 
below additional evidences of this assertion are presented. 
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Figure 6: Distributions of two principal angles (9i + fe)/2. 
Observe how curves approach the origin and touch it at 
L — 60, that indicates the two-dimensional tangencies be- 
tween contracting and expanding subspaces. The parameters 
are as in Fig. [5] 




B. Lyapunov exponents against the length of the 
system 

Figure[5]represents the Lyapunov exponents A^ as func- 
tions of L. The plots are obtained at N — 51 points of the 
spatial mesh, Ax = L/(N-1) and At w 0.01. The inter- 
val between re-normalizations and orthogonalizations is 
T/30 (see Appendix for details). Notice that the zero ex- 
ponent is absent. This is natural for the non-autonomous 
system we deal with. 

The largest exponent Ao remains almost constant as L 
varies, see the lower panel in Fig. [8] The approximat- 
ing line, obtained via least squares fit, does not have a 
noticeable slope (the slope is of the order 10 -5 ) and is 
plotted at constant value 0.138. This is equal with a re- 
markable accuracy to the theoretically predicted largest 
Lyapunov exponent ([5]) of the corresponding ODE sys- 
tem ([1]). When L is small, the system has the single posi- 
tive exponent that corresponds to spatially homogeneous 
chaotic oscillations. As L grows, the second exponent Ai 
becomes positive at L — L c . This indicates the transition 
to a spatially inhomogeneous solution. Further increase 
of L results in a cascade of passing through zero of the 
exponents. 

Fig. UK a) shows lengthes L n where corresponding Lya- 
punov exponents A n vanish. Two lines correspond to two 
sets of parameters of the system. One can see that L n 
depends linearly on n. It means that the number of pos- 
itive exponents also linearly, on average, grows with L. 
In Fig. [D(b) the intervals AL n — L n — L n _i are plotted 
(Aii = L c ). Notice that AL2 ~ ALi and these two 
values are lager then the others AL n . We attribute this 
to the transition to a non-hyperbolic attractor that takes 
place at L2. 



b) 
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Figure 7: Distributions of 6\ and (61 + 62)/%, panels (a) and 
(b), respectively, at A = 8, T — 2, e = 0.05. There are two 
positive Lyapunov exponents at L — 10, three at L — 11, and 
sixteen at L — 60. Observe the violation of hyperbolicity at 
L — 11 in the panel (a), and the emergence of two-dimensional 
tangencies at L = 60 in the panel (b). 



C. Kaplan- Yorke dimension and Kolmogorov-Sinai 
entropy 

Figure [TU] illustrates the Kaplan- Yorke or Lyapunov 
dimension -Dky [HL [Hj], and the sum of positive Lya- 
punov exponents h^, which is an upper estimate for the 
Kolmogorov-Sinai or metric entropy [bH . [l5| . Two pan- 
els are obtained for different sets of parameters. Vertical 
dashed lines mark the point L c of transition to the spa- 
tially inhomogeneous attractor, and the point L2, where 
the third Lyapunov exponent passes zero so that the at- 
tractor becomes non- hyperbolic. 

Let us consider in more details. It is known that 
for a hyperbolic attractor is equal to its Kolmogorov- 
Sinai entropy, while for a generic chaotic attractor this is 
an upper estimate for the entropy [151 ]. Because our sys- 
tem is hyperbolic at L < L2, we can use to construct 
a function which approximates the entropy at least on 
this interval. Below L c we have — Ao, while above 
this point demonstrates a power law behavior. Thus, 
employing the least squares fit, we obtain a function, ap- 
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Figure 8: Ten largest Lyapunov exponents of the system ((2} 
against L. y4 = 3,T = 5,e = 0.05. Vertical dashed lines mark 
L c w 8.44 and L2 « 15.9 (the point where A2 =0). Lower 
panel represents Ao in a large scale. The dashed approximat- 
ing line, 4 x 10~ J Z/ + 0.138, is obtained via least squares fit. 




7 computed for different parameter sets are, perhaps, 
identical (small difference can be attributed to errors of 
computations). 

The power law approximation (|15[) agrees very well 
with the numerical curve at L < L2, and at L = L2 a 
bifurcation occurs that is associated with the loss of hy- 
perbolicity. There are two possibilities above this point. 
The first one is that the Eq. (fTS"]) still gives correct value of 
the entropy, while serves as an upper estimate. The 
other possibility is that /i M correctly represents the en- 
tropy, while the approximation p5[) becomes inappropri- 
ate. Anyway, both of these variants fit well with our con- 
clusion that the system loses the hyperbolicity at L = Li- 

Above Li the entropy grows linearly with the 
length, as well as the dimension. A number of positive 
Lyapunov exponents also demonstrates a linear growth as 
follows from the linear growth of L n in Fig.OJa). This is 
a typical phenomenon for extensive fully developed chaos 
in extended systems. In particular, the linear growth of 
Dky(L) and h^L) was reported for coupled map lat- 
tices [2(|, for Kuramoto-Sivashinsky (KS) equation 1211 
and for complex Ginzburg-Landau (CGL) equation [22|. 
Also, the linear growth of £>ky was demonstrated for 
a chaotic attractor of coupled Ginzburg-Landau equa- 
tions [23l | . It can be explained by exponential decay of 
spatial correlations. Two points with space separation 
larger than the correlation length, move independently, 
so that the system can be roughly represented by a di- 
rect product of independent subsystems [2(j. Thus, the 
additivity is observed: the growth of L merely results in 
the proportional increase of the characteristic values. 



D. Spectra of Lyapunov exponents 
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Figure 9: (a) Values of length L n where corresponding Lya- 
punov exponents A n pass zero, and (b) intervals between these 
points AL n = L n — L n -i (ALi = L c ). Solid lines on both 
panels correspond to the parameters A = 3, T = 5, e = 0.05, 
and the dashed lines represent parameters A = 8, T = 2, 
e = 0.05. 



proximating as 



Ao L ^ Z/ c , 

a{L-L c y< + \ a L>L C , 



(15) 



where a = 0.083 and 7 = 0.25 for Fig. QI|a) and 
a = 0.229 and 7 = 0.26 for Fig. \Mb). The indices 



Figure II If a) demonstrates spectra of Lyapunov expo- 
nents at different L. The first curve L = 8 corresponds to 
a spatially homogeneous case when oscillations in all spa- 
tial points are synchronized and can be described by |T]). 
There is one positive Lyapunov exponent. As one can see 
from the figure, the minor negative exponents have vary 
large absolute values. It means that only a few spatial 
modes are actually involved in the dynamics, while the 
most of modes are highly damped. When L grows, more 
Lyapunov exponents becomes positive and the remain- 
ing negative exponents approach the axis of abscissas so 
that their absolute values become smaller. In the other 
words, more spatial modes participate in the dynamics. 
The separation of modes involved and not involved in the 
observable dynamics is studied in Ref. [24|. For a dissi- 
pative chaotic system is shown to exists a splitting of the 
tangent space into physical modes, responsible for the ob- 
servable dynamics, and hyperbolically isolated from them 
highly damped non-physical modes that do not bring an 
essential information about the dynamics. 

For a fully developed spatiotemporal chaos a Lya- 
punov spectrum scaled as (Dky / '/i/i)A(n/ Dky) is known 
to tend to a limiting curve at L — > 00. In partic- 
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Figure 10: Kaplan- Yorke dimension Dry, upper estimate 
for the Kolmogorov-Sinai entropy, and its power low approx- 
imation Xm l|15l) against L. (a,) A = 3, T = 5, e = 0.05. 
(b) A = 8, T = 2, e = 0.05. Vertical dashed lines mark L c 
and Z/2. 



ular, this was reported for coupled map lattice [2fJ, 
for Kuramoto-Sivashinsky (KS) equation [2lfl and for 
complex Ginzburg-Landau (CGL) equation [22]. Fig- 
ures [TTJb) and (c) represent the verification of this prop- 
erty for the system ^ at two sets of parameters. One 
can see high correspondence of curves, obtained at dif- 
ferent L. 



V. SUMMARY AND CONCLUSION 

We considered an extended system whose local dynam- 
ics is hyperbolic and spatial coupling is introduced via 
diffusion. A numerical verification of hyperbolicity of 
the attractor of this system was performed. The test 
was based on the computation of distributions of prin- 
cipal angles between contracting and expanding tangent 
subspaces of the attractor. The analysis revealed that the 
hyperbolicity is inherent only to a low-dimensional chaos 
observed at sufficiently small lengths of the system. 

The dynamics is obviously hyperbolic when oscilla- 
tions arc homogeneous in space, because each spatial 
cell merely reproduces the oscillations of a partial ODE 
system that is known to be hyperbolic. This regime is 
characterized by a single positive Lyapunov exponent. 
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Figure 11: (a) Spectra of Lyapunov exponents, and (b) scaled 
spectra at A — 3, T = 5, e — 0.05. (c) Scaled spectra at 
A = 8, T = 2, e = 0.05. 



The hyperbolicity survives when the length gets larger, 
so that the first spatial mode allowed by boundary condi- 
tions becomes linearly unstable, and the oscillations be- 
comes inhomogeneous. This transition is accompanied by 
the emergence of the second positive Lyapunov exponent. 
Further growth of the length results in the emergence of 
the third positive Lyapunov exponent. In this point the 
violation of hyperbolicity takes place. 

Beyond the point of the hyperbolicity loss, the sys- 
tem demonstrates an extensive spatiotemporal chaos that 
is characterized by a fast decay of a spatial correlation. 
We verified several standard criteria and observed be- 
havior that is typical for many others extended chaotic 
systems. Namely, the number of positive Lyapunov ex- 
ponents, the sum of positive exponents (this value is an 
upper estimate for Kolmogorov-Sinai entropy), and the 
Kaplan- Yorke dimension grow linearly against the length 
of the system. Spectrum of the Lyapunov exponents, be- 
ing properly rescaled, tends to a limiting curve as the 
length grows. 

So, if the length of the system grows and the third 
Lyapunov exponent becomes positive, we register the vi- 
olation of hyperbolicity due to the emergence of one- 
dimensional intersections of contracting and expanding 
tangent subspaces of the attractor. If the length contin- 
ues to increase, along with one-dimensional intersections, 
we observe two-dimensional ones. This is a stronger type 
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of the hyperbolicity violation, because there is higher 
probability for the perturbation to be transferred be- 
tween contracting and expanding subspaces. We expect 
that the intersections of higher dimensions also take place 
as the length diverges. It is interesting to study the viola- 
tion of hyperbolicity in the thermodynamic limit. If the 
number of modes involved in the dynamics is infinite, the 
maximal dimension of the intersections may be infinite 
too or it can have a finite value. The first case can be 
termed as a strong violation, because the capacity of set 
of the merging vectors from contracting and expanding 
subspaces is comparable with the capacity of the whole 
set of degrees of freedom. Hence, the probability for the 
perturbation to be transferred between contracting and 
expanding subspaces is non-zero. The second case can 
be termed as a weak violation. Though the intersections 
take place, the number of merging directions per degree 
of freedom is zero. Thus, the probability of the pertur- 
bation transfer vanishes. 
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Appendix: COMPUTATION OF LYAPUNOV 
EXPONENTS, COVARIANT LYAPUNOV 
VECTORS AND ANGLES BETWEEN 
SUBSPACES 

To compute Lyapunov exponents, we apply an algo- 
rithm based on the QR decomposition. See Refs. [25ll26l. 
[27J for the details of the algorithm, and Ref. [l3[ for an 
idea of the QR decomposition. 

First of all, equations for small perturbations a(x, t) 
and b(x,t) to a trajectory a(x,t) and b(x,t) of @ are 
required: 

d t a = A cos(27rf/T)5 - 2\a\ 2 h - a 2 a* - iel + d^a, 

d t b= ~Acos(27rt/T)b - 2\b\ 2 b - b 2 b* -2ieaa + dl% 

(A.l) 

where asterisk denotes the complex conjugation. To com- 
pute M\ Lyapunov exponents, we need M\ exemplars of 
the linear equation sets (|A. 1[) . which are initialized by 
an orthogonal set of random unit vectors of the length 
AN, where N is the number of points of a numerical 
mesh. Basic system @ is also initialized and advanced 
along a trajectory for a sufficiently long time to arrive at 
the attractor. Then the basic system is solved simulta- 
neously with M\ linear equation sets during some time 
interval. The more Lyapunov exponents are required, 
the shorter interval should be taken, because minor neg- 
ative Lyapunov exponents can have vary large absolute 
values so that the corresponding solutions of linear sub- 
systems decay very fast. M\ resulting vectors are then 



considered as columns of a matrix that is decomposed 
into an orthogonal matrix Q and an upper triangular 
matrix R. (An algorithm based on the Householder ro- 
tation is used [HI-) Logarithms of M\ diagonal elements 
of the R are collected, while M\ columns of the Q are 
used to re-initialize linear systems. Then this procedure 
is repeated. Averaged logarithms of diagonal elements of 
R converge to Lyapunov exponents. 

To compute covariant Lyapunov vectors according to 
the method recently reported in Ref. [12], we must do 
the similar things. After initialization of the equations, 
we make several steps no accompanied by the QR proce- 
dure, but without storing elements of R, to obtain a good 
matrix Q no . "A good" means that each linear subspace 
<S^ o , j = 1,2,.. .AN, spanned by first j vector-columns 
of Q no , contains j-th expanding (or contracting) direc- 
tion of the tangent space at uq. Starting from no, we 
make some more steps and arrive at n\. Here we have a 
matrix Q ni with columns that determine subspaces . 
Our aim now is to define arbitrary unit vectors belonging 
to these subspaces, u J n € <Sj£ , j = 1, 2 . . . AN. In fact, we 
just need to generate a random upper triangular matrix 
C ni , whose size coincides with R, and columns are nor- 
malized by 1. j-th column of C ni contains coordinates 
of v? with respect to the basis Q ni . In the other words 

U n = QnC n , (A. 2) 

where U n = {u^, u^, . . . , «„ }. Starting from C ni , we 
perform backward iterations C n -i = R^ 1 C n accompa- 
nied by re- normalization of columns of C n . Collecting 
and averaging the negative logarithms of the norms, we 
obtain Lyapunov exponents. Under these iterations the 
vectors u J n , represented by columns of the C n , are aligned 
with the most expanding directions of subspaces S^. 
These directions are associated with corresponding Lya- 
punov exponents. Because we go back in time, the high- 
est Lyapunov exponents do not dominate this alignment. 
If the number of steps from n\ to no is sufficiently large, 
getting back at no, we obtain the matrix C no with coor- 
dinates of covariant Lyapunov vectors U no , pointing ex- 
panding and contracting directions of the tangent space 
at no- Explicit form of U no can be found from (|A.2[) . 
Computed in parallel, the Lyapunov exponents allow to 
distinguish expanding and contracting directions. 

In practice, computing the covariant Lyapunov vec- 
tors for a system of many degrees of freedom, we must 
deal with very large arrays of data. For the backward 
procedure to be performed, m = n± — uq matrices R 
should be stored. The time interval between successive 
QR decompositions should be sufficiently small to treat 
minor Lyapunov exponents and corresponding vectors ac- 
curately, while the duration of the backward procedure 
must be long because the vectors are found to converge 
sufficiently slow. As a result, an array of matrices R 
runs up to several gigabytes. We recall that on 32-bit 
platforms the physical limit of an addressable memory is 
4Gb, while the memory actually available for programs 
is even less. It means that we can not store such array in 
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memory and need to write it to a file. (Otherwise, one 
can employ a 64-bit platform with appropriate amount of 
memory, of course.) Moreover, the file must be written 
in a binary format. The usual text format is not a saving 
so that an extremely large file can be obtained. 

According to Eq. (|A.2|) . we need Q nQ to restore covari- 
ant Lyapunov vectors in the original phase space. It meas 
that an array of m matrices Q must also be stored. Hope- 
fully, this is not needed. The transformation (|A.2[) pre- 
serves angles because matrices Q n are orthogonal. Thus, 
we do not need the U n to analyze the structure of the 
tangent space. Identical information about this space 
can be extracted directly from the column-space of C n . 

To compute the C n we apply a two-pass procedure. 
First, we solve the equations and perform QR decompo- 
sitions during a sufficiently long time, saving obtained 
matrices R n to a file. Then, on the second pass, we 
generate random matrix C ni , see the details above, and 
perform the backward iterations, reading R n from the 
file from the end to the beginning. When a sufficiently 
large number of transient iterations are made, we start 
to compute angles between contracting and expanding 
subspaces of the column-space of C n until arrive at the 
beginning of the file of R n . 



The algorithm of computation of the angles between 
subspaces, so called principal angles, can be found, e.g., 
in Refs. [H], [28|. Consider a matrix C n . First of all, its 
columns must be classified as vectors associated with con- 
tracting and expanding directions of the tangent space, 
according to signs of corresponding Lyapunov exponents. 
Thus we obtain a matrix S comprising of n s covariant 
Lyapunov vectors from the contracting subspace and a 
matrix U that consists of n u vectors of the expanding 
subspace. It is naturally to assume that n s > n u . For 
both of these matrices we compute the QR factorizations 
S = Q s R Sl U = Q u Ru, and then compose the matrix M: 

M = QjQ u . (A.3) 

Cosines of the sought principal angles (i = 1, . . . ,n u ) 
are equal to the singular values of the M, that can be 
easily computed, see e.g. 

This algorithm is known to fail to accurately compute 
very small angles, and in Ref. [28| an improved version is 
suggested. But, nevertheless, we use the standard algo- 
rithm, because the extremely high accuracy is not needed 
for our purposes. 
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